datum = shape of the sphere CRS = is just a reference system area and distance measurements require a projection Projection is not always needed
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(leaflet)
library(scales)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(leaflet)
ak_regions <-read_sf("data/shapefiles/data/shapefile_demo_data/ak_regions_simp.shp")
st_crs(ak_regions)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
class(ak_regions)
## [1] "sf" "tbl_df" "tbl" "data.frame"
ak_regions_3338 <- ak_regions %>%
st_transform(crs=3338)
plot(ak_regions_3338)
summary(ak_regions_3338)
## region_id region mgmt_area geometry
## Min. : 1 Length:13 Min. :1 MULTIPOLYGON :13
## 1st Qu.: 4 Class :character 1st Qu.:2 epsg:3338 : 0
## Median : 7 Mode :character Median :3 +proj=aea ...: 0
## Mean : 7 Mean :3
## 3rd Qu.:10 3rd Qu.:4
## Max. :13 Max. :4
summary(ak_regions_3338)
## region_id region mgmt_area geometry
## Min. : 1 Length:13 Min. :1 MULTIPOLYGON :13
## 1st Qu.: 4 Class :character 1st Qu.:2 epsg:3338 : 0
## Median : 7 Mode :character Median :3 +proj=aea ...: 0
## Mean : 7 Mean :3
## 3rd Qu.:10 3rd Qu.:4
## Max. :13 Max. :4
ak_regions_3338 %>%
select(mgmt_area)
## Simple feature collection with 13 features and 1 field
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -2175328 ymin: 405653.9 xmax: 1579226 ymax: 2383770
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 13 x 2
## mgmt_area geometry
## <dbl> <MULTIPOLYGON [m]>
## 1 3 (((-1156666 420855.1, -1159837 417990.3, -1161898 416944.4, -~
## 2 4 (((571289.9 2143072, 569941.5 2142691, 569158.2 2142146, 5690~
## 3 3 (((-339688.6 973904.9, -339302 972297.3, -339229.2 971037.4, ~
## 4 3 (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 654303.1~
## 5 2 (((561012 1148301, 559393.7 1148169, 557797.7 1148492, 555907~
## 6 3 (((115112.5 983293, 113051.3 982825.9, 110801.3 983211.6, 108~
## 7 4 (((-678815.3 1819519, -677555.2 1820698, -675557.8 1821561, -~
## 8 4 (((-1030125 1281198, -1029858 1282333, -1028980 1284032, -102~
## 9 2 (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186, 36741~
## 10 4 (((-848357 1636692, -846510 1635203, -840513.7 1632225, -8358~
## 11 2 (((426007.1 1087250, 426562.5 1088591, 427711.6 1089991, 4299~
## 12 1 (((1287777 744574.1, 1290183 745970.8, 1292940 746262.7, 1296~
## 13 4 (((-375318 1473998, -373723.9 1473487, -373064.8 1473930, -37~
pop <-read.csv("data/shapefiles/data/shapefile_demo_data/alaska_population.csv",stringsAsFactors = F)
class(pop)
## [1] "data.frame"
pop_4326 <- st_as_sf(pop,
coords =c("lng","lat"),
crs = 4326,
remove =F)# pass the crs that is currently is, not the one you want
pop_3338 <- st_transform(pop_4326,crs=3338)
#pop_3338<- st_coordinates() #add coordinate,
#add a colum with the ESPG
pop_joined <- st_join(pop_3338, ak_regions_3338, join=st_within)# we want to know which of the points is within each region
head(pop_joined)
## Simple feature collection with 6 features and 8 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## year city lat lng population region_id region
## 1 2015 Adak 51.88000 -176.6581 122 1 Aleutian Islands
## 2 2015 Akhiok 56.94556 -154.1703 84 6 Kodiak
## 3 2015 Akiachak 60.90944 -161.4314 562 8 Kuskokwim
## 4 2015 Akiak 60.91222 -161.2139 399 8 Kuskokwim
## 5 2015 Akutan 54.13556 -165.7731 899 1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153 777 13 Yukon
## mgmt_area geometry
## 1 3 POINT (-1537925 472627.8)
## 2 3 POINT (-10340.71 770998.4)
## 3 4 POINT (-400885.5 1236460)
## 4 4 POINT (-389165.7 1235475)
## 5 3 POINT (-766425.7 526057.8)
## 6 4 POINT (-539724.9 1456223)
pop_region <- pop_joined %>%
as.data.frame() %>%
group_by(region) %>%
summarise(total_pop =sum(population))
head(pop_region)
## # A tibble: 6 x 2
## region total_pop
## <chr> <int>
## 1 Aleutian Islands 8840
## 2 Arctic 8419
## 3 Bristol Bay 6947
## 4 Chignik 311
## 5 Cook Inlet 408254
## 6 Copper River 2294
Now we need to get it back to the correct geometries.
pop_region_3338 <- left_join(ak_regions_3338, pop_region)
## Joining, by = "region"
# you could specify, which you would need if the colums did not have the same header
#pop_region_3338 <- left_join(ak_regions_3338, pop_region, by=c("region"="region")
#dissolve the boundaries
pop_mgmt <- pop_region_3338 %>%
group_by(mgmt_area) %>%
summarise(total_pop=sum(total_pop))
plot(pop_mgmt["total_pop"])
# to keep internal boundaries
pop_mgmt <- pop_region_3338 %>%
group_by(mgmt_area) %>%
summarise(total_pop=sum(total_pop), do_union=F)
plot(pop_mgmt["total_pop"])
pop_mgmt$ESPG = "3338"
?sf::tidyverse
## starting httpd help server ... done
ggplot()+
geom_sf(data = pop_region_3338, mapping = aes(fill=total_pop))+
#geom_sf(data = rivers_3338, aes(size = StrOrder), color = "black") +
geom_sf(data = pop_3338, aes(), size = .5) +
scale_size(range = c(0.01, 0.2), guide = F) +
theme_bw() +
labs(fill = "Total Population") +
scale_fill_continuous(low = "khaki", high = "firebrick", labels = comma)
# Adding base maps in ggmap
3857 pseudomercador is the project on the tiels
pop_3857 <- pop_3338 %>%
st_transform(crs =3857)
# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")),
c("ymin", "xmin", "ymax", "xmax"))
# Coonvert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}
bbox <- c(-170, 52, -130, 64)
ak_map <- get_stamenmap(bbox,zoom = 4)
## Source : http://tile.stamen.com/terrain/4/0/4.png
## Source : http://tile.stamen.com/terrain/4/1/4.png
## Source : http://tile.stamen.com/terrain/4/2/4.png
## Source : http://tile.stamen.com/terrain/4/0/5.png
## Source : http://tile.stamen.com/terrain/4/1/5.png
## Source : http://tile.stamen.com/terrain/4/2/5.png
ak_map_3857 <- ggmap_bbox_to_3857(ak_map)
ggmap(ak_map_3857)+
geom_sf(data=pop_3857, aes(color=population),inherit.aes=F)+
scale_color_continuous(low="khaki", high ="firebrick", label=comma)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
leaflet will project your files for you but you need to make sure you use somehting lik wgs 84 (4326) in needs to be a crs with lat long
# this would be a good function
epsg3338 <- leaflet::leafletCRS(
crsClass = "L.Proj.CRS",
code = "EPSG:3338",
proj4def = "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
resolutions = 2^(16:7))
pop_region_4326 <- pop_region_3338 %>%
st_transform(crs=4326)
leaflet(options=leafletOptions(crs= epsg3338)) %>%
addPolygons(data=pop_region_4326,
fillColor = "gray",
weight =1)
#add legend scale
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),# tells how to map it to the pallette created above
weight = 1,
color = "black",
fillOpacity = 1,
label = ~region) %>%
addLegend(position = "bottomleft", #adds a legend to the map
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1) %>%
addCircleMarkers(data = pop_4326,
lat = ~lat,
lng = ~lng,
radius = ~log(population/500), # arbitrary scaling
fillColor = "gray",
fillOpacity = 1,
weight = 0.25,
color = "black",
label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m